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I \ We review our recent studies on the dynamical correlations in MC simulations 

' from the view point of the statistical dependence. Attentions are paid to the re- 

\^ \ duction of the statistical degrees of freedom for correlated data. Possible biases 

on several cumulants, such as the susceptibility and the Binder number due to 
finite MC length are discussed. A new method for calculating the equilibrium 
relaxation time from the analysis of the statistical dependence is presented. We 
apply it to the critical dynamics of the Ising model to estimate the dynamical 
' critical exponent accurately. 

Q I 1. Introduction 



Thanks to the very high performance of the computers today, we are now able to 
calculate thermodynamic quantities highly precisely using Monte Carlo (MC) simu- 
5^ \ lations. At the same time, importance of proper estimation of the statistical errors 

has widely been recognized. Since we can perform only MC runs of finite length, sta- 
tistical errors are inevitable. Suppose we have a number of statistically independent 
data; then analysis of the statistical errors are just a simple excercise of elementary 
statistics. Difficulty arises, however, when we deal with data obtained by MC simu- 
lations; since the MC simulations are dynamical methods, the data are dynamically 
correlated and thus they are no longer statistically independent with each other. 

It is still not so difficult to estimate statistical errors when we deal with such quan- 
tities as the energy and the magnetization of classical spin systems and their higher 
moments. In fact, their unbiased estimators are just the average of the measured val- 
ues at each MC step. The simplest unambiguous way for their error analysis is based 
on several independent MC runs of the same length. When a number of independent 
averages are available, then we can directly apply the elementary method of error 



estimation. For any other quantity, however, there is no guarantee in general that 
this method works properly. We can certainly get correct errors around the estimates 
by applying the error propagation law; But the estimates themselves may be suffered 
from some biases due to the finite length of the MC run. For example, estimators 
for cumulants are usually biased, although they are consistent. We can even estimate 
the magnitude of such bias, provided we know how many independent data we have. 
Problem is that what we need is not just the number of actual measurments. we 
rather need the number of statistically independent data. In fact, data taken in a 
shorter interval than the relaxation time r are strongly correlated with each other, so 
that the effective number of available statistically independent data are less (some- 
times much less) than the number of the actual measurements. Intuitively speaking, 
the elementary interval between the adjacent independent measurements would be 
2r as was discusses by Miiller-Krumbhaar and Binder some twenty years agoS. Then, 
the number of independent data is estimated as n/2T, where n is the total number 
of measurements made in each run. Based on this consideration, effect of the biases 
on the susceptibility and the specific heat - both are second-order cumulants - due 
to finite MC length was discussed by Ferrenberg et al.i Quite recently, we showed, in 
somewhat different context, that the independent data are more accurately counted 
by introducing a new scale of time, the statistical dependence time, Tdep.&S In contrast 
with the conventional reasoning where the number of independent data are simply 
proportional to the number of measurements, we found that they relates with each 
other nontrivially. 

Statistical errors and biases are not simply undesirable. Rather, they offer useful 
informations on dynamical correlations, because the dynamical correlation is just 
a different aspect of the statistical dependence. Then we can formulate yet another 
approach for studying dynamics using knowlege of the statistical errors and the biases; 
we may call it statistical dependence analysis approach. Recently, we proposed a new 
method for calculating r from a ratio of equilibrium averages of the susceptibility and 
the statistical error.i This method no longer requires calculations of time-displaced 
correlation functions; Thus its largest advantage over the conventional methods is 
that it enables us to make unambiguous statistical analysis, because we can simply 
follow the standard methods used for static quantities. We have applied this method 
to the MC simulations of the two- and three-dimensional Ising model, and calculated 
the dynamical critical exponent z highly accurately. 

2. The Statistical Dependence and dynamical correlation 

First, we discuss the relation among the susceptibility (second-order cumulants, in 
general), the relaxation time and the statistical error. In the course of the discussion, 
we will examine the method for counting the number of the independent data properly. 

Suppose we make N identical MC simulations which differ only in the random 



number sequences. Random number sequences should be independent from run to 
run, in order mutual statistical independence of these simulations to be guaranteed. 
We make n measurements of a quantity Q in each run. Let us define two types 
of averages as follows: (l)The average over the data taken in a single run, {Q)a = 
^ Ylk=i Qa{k), and (2)The average over all the independent runs, (Q) = Y.a=i{Q)a, 
where Qa{k) denotes the value of the quantity Q at the k-th measurement in the a-th 
run, with k = 1, 2, ■ ■ ■, n and a = 1, 2, ■ ■ -, A^. 

Elementary statistics says that when all the n data in a run are statistically inde- 
pendent, then the statistical degrees of freedom (DOF) also is n and the expectation 
value for the variance of (Q) (we omitted the suffix a here by an apparent reason) is 
simply given as 

E{((Q)-/i)^} = ^, (1) 

where /i = E {Q} and = E {Q^} — [J? are the expectation values of Q and of 
the variance of Q, respectively. In the present situation, however, the data in a run 
are dynamically correlated with each other rather than they are independent. The 
reduction in the DOF due to these correlations can be discussed in terms of the time- 
displaced correlation function of Q. The expectation value for the variance of {Q) 
can be written as 

E {((g) - lif] = - E (e {m{k)m{l)} - = ^ E E {C{k, I)} . (2) 
^ k,i=i ^ k,i=i 

where C{k, I) is the (normalized) time-displaced correlation function of the fluctuation 
in Q for "time" k and /. By comparing Eq. (^ with Eq. (p, we find that the DOF 
is reduced from n to rired = I YJki=i E {C{k, I)}. Thus the reduction factor n^ed/n 
depends highly nontrivially on both n and r. 

Let us introduce a new quantity, rjep with the dimension of time as 

1 " 

^dep^^ E E{C(fc,/)}. (3) 
k,l=l 

Then we can put Eq. (H) in a similar form as Eq. ([l|): 

E{((Q)-/i)^} = ^a^ (4) 

The DOF is thus reduced by a factor l/2rdep due to the dynamical correlations. In 
other words, we can interprete 2rdep as the mean interval of successive statistically 
independent measurements. It may be somewhat surprising result, because this in- 
terval depends not only on the relaxation time r but also on the total number of the 



measurements n; Namely, this interval is determined only after all the measurements 
have been made, even if we know r beforehand. 

Equation (^) tells us how we can estimate r^ep and, as a result, the reduced DOF 
by simulations: The unbiased estimator {6Qy for E {{{Q) — /i)^} is calculated from N 
independent runs as {6Qy = ((Q)^ — {Q) ) with \6Q\/\/N being the statistical 
error in Q. The variance cr^, on the other hand, is calculated as the susceptibility xq 

2 

associated to Q, xq = (Q^) ~ {Q) ■ Thus we get the estimator of r^ep as follows: 



_ n(SQr- _ nN {Qy - (Q) 
2X0 -2(iV-l)(g^_(Qy- 

where the r.h.s. can be calculated by simulations. Let us see how Tdcp behaves. For 
simplicity, we assume a single exponential relaxation for the correlation function, that 
is, we put E{C{k,l)} = exp(— — /|/r). If we approximate the summation in Eq. (H) 
by an integration in the range [0, oo), we reproduce the result obtained by Miiller- 
Krumbhaar and BindeiS, which implies that the reduced DOF is given as rired = n/2T. 
From the derivation, however, this expression is valid only if two conditions n <^ t 
and r <^ 1 are satisfied, that is, in the long-time and the long- relaxation-time limit. 

We can go a little further. The summation in Eq. (0) can be calculated explicitly 
without assuming the long-time nor the long- relaxation-time limit. And we get 



''"dep 



1 + A 2A(1-A" 



1-A 72(1 -A)^ 



(6) 



where A = e"^''^. It can easily be confirmed that the reduced DOF, rired = njljAep-, 
varies smoothly from 1 for n = 1 to n/2r for n — oo as it should do. Moreover, we 
can use Eq. (||) for estimating r, once Tdep is estimated by the simulations. 



3. On Unbiased Estimators of Cumulants 



As we have mentioned in the introduction, estimators for cumulants of any order 
are, in principle, biased. For example, it is well known that the unbiased estimator for 
the variance, that is, the second-order cumulant, from n independent measurements 
is given as 

XQ = ^ {{Qy - iQr) ■ (7) 

The factor n/{n — 1) appears here for correction of the bias due to the finite number 
of the measurements. Therefore, the estimators for susceptibilities calculated from 
the fluctuations in general should also be corrected by the above factor. Otherwise, 
the calculated susceptibility will be systematically underestimated. Importance of 
correcting these systematic errors on susceptibilities was discussed in ref. 2. In actual 



Fig. 1. n dependence of the magnetic susceptibility for 16"^ 3D Ising model. The solid line indicate 
Eq. (I). 



simulations, the measurements are not independent of each other, as has been dis- 
cussed so far. Therefore, we can not simply take n/{n — 1) as the correction factor. 
In ref. 2, it is argued that the susceptibility estimated from n measurements Xq(^) 
and its true expectation value xq(oo) is related as 

XQin) = xq(oo) (l - . (8) 

This form coincides with the actual MC data for rather large n. It is not surprising 
that eq. @ deviates from the actual MC data for smaller n, because n/ (2r + 1) was 
used as the DOF, which, as was discussed in the previous section, is valid only in the 
long-time limit. 

As we have shown, the proper form of DOF is n/2T^ep rather than ri/2r. Therefore, 
Eq.(||) should be modified as 

XqH=xq(oo)(i-^). (9) 

For examining this equation, we made MC simulations of three-dimensional Ising 
model at the criticality taking several different lengths of the run. The system size is 
16^. Figure |I] shows n dependence of the magnetic susceptibility x calculated as the 
second-order cumulant. The solid line is the expected behavior from Eq. (1); xq(cx3) 
and T were estimated from a longer run, and then rjep was estimated using Eq. (||). 
We can clearly see the systematic deviation of the shourt-run values of x from the 
long-time average. Moreover, this systematic behavior is really expressed by Eq. (^) 
even for rather short simulations. By making the similar plot for the specific heat C, 
we found that that C also behaves as expected. Thus the validity of eq. (||) has been 
verified. 

So far, we have examined the bias on the second-order cumulants. The origin of 
the bias is attributed to the fact that the second-order cumulants are calculated as a 
combination of the estimators of the first-order and the second-order moments, and 
thus are calculated only after the averages of all the involved moments are taken. 



Fig. 2. n dependence of (m'')/(m^)^ for 16'^ 3D Ising model. 



Another important quantity frequently used in study of phase transitions is also in 
this category: that is, the Binder number. The relating cumulant to the Binder 
number is a simplified version of the fourth-order cumulant, 3{Q'^)^ — (Q^). It is 
easily confirmed that its unbiased estimator is given as follows: 

3E{{Qr} - E {(Q^)} = 3^(a^ + + ^/Z4, (10) 
where = E{Q^}. By dividing it by E{((5^)^} we get 

E{(g^)} _ 

E{(Q2)n l + i(i?oo-l)' 
where i?oo is the true value for the ratio of the moments: 



_ E{gn _ .... 



It, however, is not what we can use for deriving the unbiased estmator for the Binder 
parameter; while it is the ratio of the expectation values of the moments, what we 
really need is the the expectation value of the ratio, E • Only if the 

correlation between the denominator and the numerator can be ignored, we can use 
Eq. (|TI|) as its approximant. That is certainly not true in general, of course. But, in 
any case, the bias of 0{l/n) is expected for the Binder number. 

In Fig. ^ n dependence of the ratio {m^) / {rn?)'^ for the magnetization moments 
is plotted. The effect of the bias is clearly seen. But Eq. (lTl|) cannot reproduce this 
behavior. The similar plot for the energy, {E^)/{E'^)^ is shown in Fig. ^. The solid 
line indicate Eq. ([Tl|) . In contrast with the case of the magnetization, we see that 
the behavior of the ratio is well approximated by Eq. (|ll|)- This difference between 
two Bender numbers may be attributed to the fact that the Binder number was 
originally derived from a simplified version of the fourth-order cumulant, with the 
first- and the third-order moments set to be zero beforehand. For the magnetization, 
this simplification is meaningful, since the expectation values for all the odd-order 



Fig. 3. n dependence of (E^) / (E^)^ for 16'^ 3D Ising model. The solid line indicate Eq. (11). 



moments vanish due to the time-reversal symmetry. The energy, on the contrary, 
does not have this property. Therefore, the energy Binder number has only a vague 
relation with the fourth-order cumulant. That may be the reason why the behavior 
of the energy Binder number agrees with Eq. (]TT|), while the magnetization Binder 
number does not. In any case, their estimators approaches the true values slowly with 
n. 

For the quantities we have dealt with, that is, the second-order cumulants and 
the Binder number, we can easily realize that their estimator should be biased. Sup- 
pose we have only one measurement, namely, take n — 1; their estimators give just 
trivial values in this extreme limit: for the second-order cumulants and 2/3 for 
the Binder number irrespective of the measured values of the magnetization or the 
energy. Therefore, it is obvious that the values of these estimators vary from these 
trivial values to the true values as n increases. 

4. Estimation of the Dynamical Critical Exponent 

As we have discussed, we can estimate the relaxation time r using the averages of 
the statistical error and the susceptibility. Since only imformations of static averages 
are used in this calculation, it allows us of an unambiguous error estimation for r. As 
a result, we can obtain a highly reliable estimate of r. Thus it is a suitable method 
especially for calculating the dynamical critical exponent z. We have applied this 
method to two- and three-dimensional Ising models at their criticality, and obtained 
very accurate estimates of z. In what follows, we present the results briefly. Details 
of the calculations are found in ref. 3. 

Figure ^ shows the log-log plot of the relaxation time r against the linear dimension 
L of the three-dimensional Ising model. We estimated r from magnetic susceptibility 
and the statistical error of the magnetization using Eq. (^). The data fit very well to 
a straight line. According to the dynamical finite-size scaling theory, the slope of the 
plot gives z. After careful statistical analyses, we have got z = 2.03 ± 0.01. Seeing 
the statistical error, the value of z we obtained is the most accurate one proposed so 



Fig. 4. The finite-size scaling plot of r for 3D Ising model. The slope of the line gives z = 2.03±0.01. 

far, as far as we know. By making the similar finite-size scaling analysis, we have got 
z = 2.173 ± 0.016 for two-dimentional Ising model. The present statistical error is 
also one of the smallest ones among the studies so far. Most importantly, we did not 
use any special technique for estimating the statistical error in z] We just followed 
the standard error analysis methods used for static quantities. 
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